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Abstract 

In the paper, a novel algorithm employing pseudo-spectral approach is developed for the PKN mo- 
del of hydrofracturing. The respective solvers compute both the solution and its temporal derivative. 
In comparison with conventional solvers, they demonstrate significant cost effectiveness in terms of 
balance between the accuracy of computations and densities of the temporal and spatial meshes. 
Various fluid flow regimes are considered. 

1 Introduction and preliminary results 

Hydraulic fracturing is a widely used method for stimulation of hydrocarbons reservoirs. This technology 
has been known and successfully applied for a few decades [HI [T2j |6]. Recently it has been revived, 
due to economical reasons, as a basic technique for exploitation of unconventional deposits of oil and 
gas. The phenomenon of a fluid driven fracture propagating in a brittle medium is also present in many 
natural processes (e.g. magma driven dykes - [31], subglacial drainage of water - [34]). 

Throughout the years, starting from the pioneering works of Sneddon and Elliot [35], Khristianovic 
and Zheltov [T3] , Perkins and Kern [3D] , Geertsma and de Klerk [TT] , and Nordgren [5^ , various models 
of hydrofracturing have been formulated and used in applications. A broad review of the topic can be 
found in |15[|16[[^ll7j . Together with increasing complexity of the models describing this multhiphysics 
process, the computational techniques have been continuously enhanced. A comprehensive survey on 
the algorithms and numerical methods used in hydrofracturing simulation can be found in [4] [9] . 

Responding to the recent demand, an increasing stream of publications have appeared concerning 
additional information on seismic events, shear stresses in the rock formation, multifracturing and others 
and their implementation into the solvers (35] [521 HT] |H] ■ Also, a considerable effort has been made to 
improve the existing algorithms by incorporating new efficient numerical techniques [18] [21] [28] [29] . 

The main computational challenges associated with the modelling of hydraulic fractures are: a) 
strong nonlinearity resulting from the coupling between the solid and fluid phases, b) singularity of the 
gradients of the physical fields near the crack tip, c) moving boundaries, d) degeneration of the governing 
equations at the crack tip, multiscaling and others. 

To achieve the maximal possible efficiency of numerical simulations, the computational algorithms 
should be formulated in proper variables accounting for all the problem peculiarities [21j . As a result, 
they allow one to reduce the volume of processed data, which is especially important when dealing with 
complex geometries and/or multifracturing. 
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The analysis presented in this paper is devoted to the PKN model of hydrofracturing. This model 
contains all the peculiarities mentioned above, except for the non-local relation for the fluid-solid cou- 
pling. Although we restrict our interest only to a single fracture, the developed algorithms, thanks to 
their robustness, can be successfully applied to model a system of cracks. 

The numerical analysis of the problem should be backdated to Nordgren [55] who extended the 
Perkins and Kern model [50] to account for the fluid loss effect and fracture volume change. As a result, 
the crack length was determined as part of the solution. The author proposed a finite difference scheme 
to solve the problem, which is in fact equivalent to the finite volume (FV) method. 

Further development of the PKN formulation was done by Kemp fl3| . who (a) implemented the 
specific boundary condition at the moving crack tip into the FV scheme, (b) incorporated asymptotic 
behaviour of the solution near the crack tip in a special tip element, (c) indirectly used the fourth power 
of the crack opening (w^) as a new dependent variable, instead of the crack opening itself. For the 
early-time asymptotic model Kemp proposed a power series solution, presenting its first four terms. 

The recent paper by Kovalyshen and Detournav|15j has extended most of Kemp's results, incorpo- 
rating all information on the PKN model available to date. They present various asymptotics, complete 
analytical solution for an impermeable rock (directly extending the results from |13| from four leading 
terms to an infinite series representation), FV algorithm with a special tip element and a numerical 
benchmark for the Carter leak-off, linking the results to the scaling approach developed in [3j [71 [23l [24] . 

In [in]-[ll], the PKN model was reformulated by Linkov to improve the efficiency and stability of 
computations by (i) introducing proper dependent variables (cubed fracture opening, w^), (ii) utilizing 
the speed equation and (iii) by imposing a modified boundary condition at a small distance behind 
the crack tip (e-regularisation) . Additionally, the analytical solution for an impermeable rock was 
evaluated for the new dependent variable in a form of rapidly converging series in |21| . Moreover, the 
author highlighted in [19] that numerical schemes exploiting a fixed position of the crack tip during the 
iterations may become ill-posed. 

In |22j and |17j the e-regularisation technique was further enhanced by (i) appropriate adaptation of 
the speed equation to the chosen numerical scheme and (ii) improved way of imposing of the regularized 
boundary condition. A detailed discussion on various aspects of application of implicit and explicit 
numerical schemes was provided. 

In this paper we are presenting a novel algorithm based on the pseudo-spectral approach. Namely, 
we propose an efficient numerical algorithm to solve a specific self-similar problem and extend the results 
to the general (transient) formulation. Since, the integration schemes used in the algorithm incorporate 
the exact boundary conditions at the crack tip, no regularization technique is necessary. The most 
accurate two points representation of the temporal derivative is used to guarantee an optimal algorithm 
performance. Finally, two solvers are developed which show their robustness and stability. They both 
demonstrate high cost effectiveness in terms of the relationship between the volume of the processed 
data and the accuracy of computations. Moreover, additionally to the crack opening and length, the 
temporal derivative of the former and the crack tip velocity are automatically returned as components 
of the problem solution. 

1.1 Problem formulation 

Let us consider a symmetrical crack of length 21 situated in the plane x € [—1, 1]. The crack is fully filled 
by a Newtonian liquid injected at the middle point {x = 0) with a known rate qo{t). Note here, that the 
crack length evolution, I = l{t), is the result of fiuid fiow inside the fracture. Due to the symmetry of 
the problem, one can restrict the analysis to the half of the crack x G [0, l{t)]. 

The classic mathematical formulation of the PKN model of hydrofracturing was given in [26] . Below 
we present a system of equations constituting the model. The mass conservation principle is expressed 
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by the continuity equation: 



^ + ^ + g/ = 0, t>to, 0<x<l{t), (1) 



while the Poiseuille equation describes the flow in a narrow channel. In the case of a Newtonian fluid, 
it is written in the following form: 

Here w = w{t, x) stands for the crack opening, q = q{t, x) is the fluid flow rate, p = x) {p — pf ~ ao, 
Co - confining stress) refers to the net fluid pressure. The constant M , involved in the Poiseuille equation, 
is computed as A/ = 12/i, where ^ denotes the dynamic viscosity (see for example [1]). The function 
qi ~ qi {t, x) from ([1]) is the volumetric rate of fluid loss to formation in the direction perpendicular to the 
crack surfaces per unit length of the fracture. This function is usually assumed to be given, but it may 
depend on the solution itself as well. To account for various leak-off regimes, we accept the following 
behaviour of qi : 

qi{t,x) = Qi{t){l{t)-x)\ for x^l{t), (3) 

for some constant -q > —1/2. Note that the case ?/ = —1/2 corresponds to the Carter law [5], while 
^7 ^ 1/3 guarantees that the leak-off vanishes near the crack tip as fast as the crack opening at least (see 
for details [17]). 

The group of fluid equations is to be supplemented by the relation describing deformation of the rock 
under applied hydraulic pressure. In the case of the PKN model, a linear relationship between the net 
fluid pressure and crack opening is in use: 

p = kw, (4) 

where a known proportionality coefficient k = -^j^^ is found from the solution of a plane strain 
elasticity problem [26] for an elliptical crack of height h. E and v are the elasticity modulus and 
Poisson's ratio, respectively. 

The above equations are equipped with the boundary condition at a crack mouth {x — 0) determining 
the injection flux rate: 



k 

~ M 

and two boundary conditions at a crack tip: 



3 ^""^ 
dx 



qo{t), (5) 



x=0 



w{t,l{t))^0, q{t,l{t))^0. (6) 

In order to define the crack length, l{t), the global fluid balance equation is usually utilized (see for 
example [1]) 

ri(t) rt 

/ [wit, x) — w{0, x)]dx — / qo{t)dt+ 

Jo Jo /-tN 

qi{t, x)dtdx = 0. 



Finally, the initial conditions are assumed in the following way: 

w(0,x)=0, ?(0)=0. (8) 

System (d} - (|5]) constitutes the classic formulation of the PKN problem. It was shown in [T3] and 
|10] that the asymptotic behaviours of w and q near the crack tip are interrelated, and the first term of 
the expansion for the crack opening may be written as: 

w{t,x) wo{t){l{t) - x)°' , as x-)-l(t). (9) 
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For the classic PKN model the exponent a — l/'i was found in [T3]. Thus condition ©2 is always 
satisfied as it follows from ([2]) and ([4]). As a result, the model does not account for the standard stress 
singularity of fracture mechanics at the crack tip, and thus is relevant for the so-called zero toughness 
regime (see e.g.[T]). 

Remark 1. Despite that zero crack opening and length arc considered as the initial conditions, 
all authors begin their studies from the asymptotic model for the small time. With the assumption of 
zero leak-off term in the continuity equation and constant go, the problem is reduced to a self-similar 
formulation. The full numerical analysis is then continued by taking the similarity solution as the initial 
state. This effectively means that the initial conditions (|8]) can be replaced by the non-zero crack opening 

= /o, u)(0, x) = w^{x), X £ (0, lo). (10) 

In this paper, the modified formulation of the PKN model is considered, following the recent advance 
in the area of numerical modelling [121 HOI HIl HI] • Thus, to trace the fracture front we use the so-called 
speed equation, instead of the fiuid balance relationship ([T]): 

^ = ^'o(t) = |L.,,. (11) 

The speed equation assumes that the fracture tip coincides with the fluid front, which excludes the 
presence of a lag or an invasive zone ahead of the fracture tip. Originally it was introduced by Kemp 
[T3] and has been recently revisited by Linkov [191 ISHl HI] . 

Note that, on substitution of equations ([2]), ([4]) and ([9]) into (fTT|) . one obtains a relationship between 
the crack propagation speed and the multiplier of the leading term of the crack opening asymptotic 
expansion 

i = ^fe"'"*"- <'^) 

This implies that the quality of the numerical estimation of Wq ( see estimate ) should he vital for the 
accuracy of computations. 

By substituting the Poiseulle equation ([2]) into the continuity equation (|T]) one obtains a lubrica- 
tion (Reynolds) equation for the considered problem, where the net fluid pressure function p{t, x) is 
eliminated: 

In this way the modified formulation of the PKN model includes: i) the Reynolds equation p^ : ii) 
the boundary conditions (0) - (Oi; iii) the asymptotics (jH]); iv) the initial conditions (|10p : v) the speed 
equation in the form ([12]). 

The paper is organized as follows: in the next subsection we present the normalized formulation 
of the problem. Then, two types of self-similar solutions for the PKN model are discussed. These 
solutions are used in section 2 to investigate a numerical algorithm for a time independent variant of the 
problem. In section 3, the algorithm is modified to tackle the transient regime. Two alternative integral 
solvers are developed and their performances and applicability are examined. Section 4 contains the 
final conclusions. 

1.2 Normalized formulation. 

Following |17j . we normalize the problem by introducing dimensionless variables: 

X ~ t M , . 
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,~ , wit, x) l(t) o , , , , 

w{t,i)^^-^, L{t) = ^, Z^goW-MoW, (14) 



where x S [0,1], L(0) = 1. 

In the new variables equation ((T3)) reads: 



dw _ L' dw 

X- 



dt L dx L^{t) dx\ dx 

f > 0, < i < 1. 
The boundary conditions ([5]) - (|6l)i may be rewritten as: 

1 



(15, 



~3 

dx 



^qoit), w{t,l) = 0. (16) 

x=0 



m 

The initial conditions pO|) are defined as: 

L(0) = 1, wiO,x) = w^i), ie[0, 1]. (17) 
The asymptotic expansion for crack opening Q takes the form: 

i&(i,i) ~^&o(^)(l-^)^^^ for i ^ 1. (18) 

For the sake of completeness of the normalization, we also present the global fluid balance equation 
([7]), although it will not be used later on: 

L{i) I w{t,x)dx— / w{0,x)dx— / qi){t)dt 

Jo Jo (19) 



+ / L{t) / qi{t,x)dxdt = 0. 

^0 "'0 

Finally, the transformation of the speed equation yields: 

As shown in [52], equation ([50)) is convenient to trace the fracture front when standard ODE solvers 
are in use for the dynamic system describing the problem. On the other hand, the crack length can be 
computed from ([501 by direct integration to give: 
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L(i) = ^l + ^J^wUr)dT, (21) 

which is useful when an implicit method (for example Crank- Nicolson scheme) is utilised (see |22|). 
On substitution of (l20l) into dTSl) one can rewrite the later to obtain: 



3L^(^ + ,.(^,^))--o^^+3-(.3_j. (22) 

In the following, equation ([22]) will be used as a basic relation to formulate our integral solvers. 

From now on, for convenience, we shall omit the tilde symbol in all quantities. In this way all the 
notations refer henceforth to the normalized formulation. 
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1.3 Self-similar solutions 



Let us assume 

qiit,x)=je'y'q:ix), (23) 
and look for the similarity solution of the problem in the form; 

w{t, x) = u{x)e'^\ wo{t) = uoe^\ (24) 

where the asymptotic behaviour (|18[) holds true, and mq is the limiting value of u defined in the same 
manner as in estimate (|18p . Thus, equation (|20p transforms to an identity if one takes: 

L^(t) = ^e3^*. (25) 



On substitution of (|23p , \2A\ , and (p5|) into the equation ((22)) one can reduce the latter to the following 
ordinary differential equation: 

Pul{u + q*i)=A{u), (26) 

with (3 = 2/3. Here, the nonlinear differential operator A is defined by the right-hand side of equation 
and is equipped with the boundary conditions 



Q -3/2 



3 diL 
dx 



= q*0, u(l) = 0, (27) 

x=0 



where we have introduced an auxiliary notation: 



9o = y|e-^go(t). (28) 

If ^0 is constant, then equation (|26p together with the boundary conditions (P7)) do not depend on time 
and constitute a boundary value problem (EVP) degenerated at point x = 1. Indeed, the nonlinear 
coefficient in front of the second order term of the differential operator vanishes at the point x = 1 in 
view of the boundary condition (P7)) 9. This EVP is in fact a self-similar formulation of the original 
problem with specific, given leak-off regime and the inlet flux. 

Other class of similarity solutions can be found, for some a > 0, by assuming: 

qi{t, x) = ^{t + a)'^^^q*{x), w(t,x) ^ {t + a)^u{x), 

wo{t) = uo{t + ay , (29) 

As a result, one again obtains the EVP ^ - ^ with /3 = 27/(87 + 1) and 



'?o* = ^(^(« + t)'-^'Zo(0- (31) 

Thus, for 7 = 1/5 the self-similar solution corresponds to the constant injection flux rate, while the crack 
propagation speed decreases with time as L'{t) = 0{t~^^^) for t ^ 00. If, however, one takes 7 = 1/3, 
the crack propagation speed is constant and the injection flux rate increases with time: qo{t) — 0{t^/'^) 
for t — !■ 00. 



6 



Note that self-similar solutions do not necessarily satisfy the initial conditions P7|) as the normalised 
initial crack lengths are: 

for the first and the second type, respectively. 

Remark 2. As one can see the second type of similarity solution has a physical sense for any 
— 1/3 < 7 < 00 and thus, can be used to model three different transient regimes of the crack evolution: 
crack acceleration (7 > 1/3), crack deceleration (7 < 1/3), and a steady-state propagation of the fracture 
(7 = 1/3). The first type of solution possesses a physical interpretation only for positive values of 7, 
which restricts its application to the cases of accelerating crack. 

The self-similar solutions formulated above are used in the following sections to analyse computational 
accuracy provided by the developed solvers. So far to this end, the asymptotic models have been usually 
employed [26l [131 [151 [21] . However, all of them are restricted to the case of a constant influx, go- 



2 Numerical solution of the self-similar problem 

In this section we will formulate an algorithm of the solution for the self-similar problem defined by 
equation ([26)) and the boundary conditions ([27)) . The following representation of the sought function 
u{x) will be accepted: 

u{x) =uoil-xy^^ + Au{x). (32) 

It results from the asymptotic behaviour ([TS)) and Au{x) ~ 0((1 — x)'^) for x — ?► 0. Parameter ^ > 1/3 
depends strongly on the behavior of the leak-off function qi near the crack tip. In particular, when qi 
vanishes near the crack tip in the same manner as the solution, or faster, {rj > 1/3) then C = 4/3. One 
can show that (compare ([3])) 

C = min{4/3,l + ry} > 1/2, (33) 

see also [T7] for details. 



2.1 Integral solver for the self-similar problem. 



Below we present an algorithm to solve equation ([26| by numerical inversion of the operator A. Exploiting 
the solution representation ()32p . the inverse operator defines both components: uq and Ait. To 
derive A~^, we integrate the equation ()26p twice over the interval [x, 1] taking into account the boundary 
condition ^7^ 9. Then, after simple transformations one obtains: 



3ug(l - x)Au = -|[6u§(l - x)2/3(Am)2+ 
Au„{l-xy/^{Au)^ + {Au)^]+ul f Aud^+ 



2ul I (^ - x)ud^ - (1 - x)ul I ud^+ ^^^^ 



J X 

In short, the latter can be symbolically written in the compact form: 



Au = Gi (13, uo, Au) + Wo, Aw, 9/*), (35) 
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where the operators involved in the right-hand side are defined as foUows: 

3^3(1 - x)Gi = -| [6(^2(1 - x)2/3(Aw)2 + 

J X 



(36) 



3(1 - x)G2 I A 

J X 



udi + -{?,p-\)uo(\~xy'^+ 



/3 / (e - x)q*idi. 

J X 

One can conclude from ([3]) and (fT5)) that for x — )■ 1 

G,=0{{l~xY+^),G2 = 0{{l^x)^) 



(37) 



(38) 



where C, defined in ((33)) . 

The relation to compute mq is derived by integration of (|26p with respect to x from to 1. Then, 
taking into account the representation (15^ and the boundary condition (P7|) i . one can formulate the 
condition: 

G:i{P,uo,Au,q*)=Q, (39) 

where 

5/2 



3/2 



G3 = -(/3 + lX^+ 



(/3 + 1) / Audx + /3 / qidx 



%■ 



It is easy to prove that for any gg > and /? > — 1 there exists a unique positive solution uo of equation 
([39]) . regardless of the values of the functions Am and gf. 

The inverse operator is defined, by equations ([55]) and ([5^ . as: 



[uo,Au]^A-\l3,ql,q*). 
Its numerical execution is based on the following iterative algorithm: 

0, 



(40) 



'G3 (/3,4^+^\AMW,gS 

G2 (/3,4^+^\AuW,gr) 



(41) 



where superscripts refer to the consequent iterations. In the first step we assume that: 

Au^°^ = 0. 



(42) 



Of course, if any better approximation is available, it can replace (P^ . Note the first relation of (|^T|) is a 
nonlinear algebraic equation that can be solved e.g. by the Newton-Raphson method, while the second 
equation is a typical iterative relationship. 
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2.2 Numerical examples and discussions 

In this section we investigate the performance of the numerical algorithm (|4ip - (|1H). To this end, 
the power law type self-similar solution (|29p is utilized as a benchmark. Apart from the fact that the 
equations for other self-similar solution (j24[) look identically, the value of the parameter f3 appearing 
in the exponential benchmark is always the same (/3 — 2/3). Thus the power law type self-similar 
benchmark gives us an opportunity to manipulate with the value of /3, which will be crucial for further 
implementation of the algorithm to the transient regime. 

Let us utilize the following variant of benchmark solution used previously in |22| (compare ()63|) in 
the Appendix A): 

= + (43) 

.(x) = -^Q-^) (1 -X) +0.05(1 -.f, 

which yields 77 = 4/3. 

Computations are carried out for different meshes based on iV + 1 nodal points: 

=1- (1-^)', j = 0,l,...,iV. (44) 

When one sets g = 1, the points of spatial mesh are uniformly distributed over the whole interval - 
this mesh will be called the uniform mesh. By taking g > I we obtain increased mesh density while 
approaching the end point x ~ I (the larger g the greater mesh density near the crack tip) - this mesh 
will be refereed to as the non-uniform mesh. In our computations wc will use g = 3. This choice is 
motivated by the following reason. To increase the solution accuracy, we compute integral operators 
from ()4ip using the classic Simpson quadrature, which gives an error controlled by the fourth derivative 
of the integrand. Accounting for the asymptotic behaviour of the solution, the transformation (j44p . for 
£1 = 3, improves the smoothness of the integrand with respect to the new independent variable near the 
crack tip (replaces the fraction powers function). In this way, the error of integration can be minimized. 
We shall confirm this below in numerical tests. 

We use two parameters as the measures of computational accuracy: (?) - the maximal point-wise relative 
error (denoted as: Su) of the complete solution u{x) and (ii) ~ a relative error (denoted by Suq) of the 
coefficient uq defining the leading asymptotic term in p2p . 

The efficiency of computations will be assessed by the number of iterations needed to compute the 
solution. The iterative process is stopped in each case when the L2-iiorm of the relative difference 
between two consecutive approximations becomes smaller than e = 10"^*^. Finally, we analyze how the 
solution errors depend on the number of nodal points N. 

The computations revealed that convergence of the iterative process (|^T|) may only be achieved for 
some range of (3 values. For the analyzed benchmark it was: /3 G [—1.8,4.8]. In fact, this interval is 
wider than one could expect and fully covers any physically motivated values of /3 G (0,2/3) following 
from the self-similar formulation. Interestingly, there are also solutions for /3 < — 1 (compare discussions 
after equation ([39])). This information shall be used later on, to construct a solver for transient regimes. 

In Fig. the values of 6u and Suq obtained for the analyzed benchmark ([43]) are shown. The 
computations were done for the meshes composed of 100 nodes. As can be seen, the non-uniform mesh 
gives at least one order better accuracy of the solution, for the same number of nodal points. The error 
of uq is lower than the error of the complete solution u(x), as one could expect, while its distribution 
is not as smooth as for 6u(x). The minimum of 6u(x) is located near /3 ~ 1/3 which corresponds to 
the steady-state similarity solution. Interestingly, in the case of the uniform mesh the aforementioned 
minimum is deeper and sharper than for the non-uniform one. 

Fig. [TId) depicts the number of iterations needed to obtain the final solution for different values of 
j3. It shows that the convergence rate almost does not depend on the type of mesh chosen. The best 
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Figure 1: a) The accuracy of computations: maximal relative error of u{x) and uq as a function of /3; 
b) Number of iterations to reach the final solution as a function of /3. Dashed lines refer to the uniform 
mesh, solid lines to the non-uniform mesh (for g ~ 3). 
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Figure 2: The maximal relative error of solution as a function of the number of nodal points N: a) 6u, 
b) Suq. Dashed lines refer to the uniform mesh, solid lines to the non- uniform mesh (for g = 3). 



efficiency of computations appears for approximately |/?| < 1/2. Note that this interval corresponds to 
the values of /3 which provide the best accuracy of computations. The time of computations (number 
of iterations) increases with growth, however this trend does not exhibit a bilateral symmetry (for 
negative and positive values of /3). Moreover, the cost of computations increases with the solution error. 

Fig. [2] shows how the accuracy of computations depends on the number of nodal points A^. Both 
chosen meshes are analyzed for three different values of /3 ~ —1, 0, 1. It can be seen that the non- uniform 
mesh gives at least two orders of magnitude better accuracy than the regular one. For both types of 
meshes, the values of Suq are much lower than respective Su, however the best result for Suq does not 
necessarily correspond to the best Su (e.g. /3 = — 1 for the regular mesh). The non- uniform mesh 
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Figure 3: The maximal relative error of solution 5uq for various spatial meshes. Computations were 
done for P — 1/3. 

provides much lower sensitivity of solution accuracy to the variation of /? than the uniform one. Also 
the maximal level of accuracy is obtained much faster for the non-uniform mesh. In the considered case, 
it is sufficient to take only 60 nodal points to achieve the maximal possible accuracy. 

The last test in this subsection identifies the influence of spatial meshing on the solution accuracy. 
To this end, we consider the following values of g = 1, 2, 5 from representation ((44|) . The benchmark 
taken here accepts /3 = 1/3, as it provides the best accuracy and will be important in next subsections. 
For each of the values of p, a characteristic 5u{N) was computed. The results are presented in Fig. [3l 
It shows that, regardless of the mesh under consideration (or equivalently, the value of the parameter 
q), the maximal achievable accuracy is the same. However, this ultimate level is reached for different 
numbers of nodal points, A^. The fastest convergence takes place for (? = 3, which confirms our previous 
predictions on the optimal choice of the spatial meshing. The slowest convergence to the saturation level 
manifests the uniform mesh [g — 1). 

The overall influence of the value of the parameter g on the accuracy of computations results from 
the following trend: the larger the value of g the lower error near the crack tip and the greater error 
near the crack inlet. The optimal balance between the local errors is observed for ^? = 3, which confirms 
our predictions. In the following, only the non-uniform mesh for g = i will be used in computations. 

3 Solution to the transient problem 

In this section we adopt the idea of the integral solver, developed for the self-similar formulation, to 
the transient regime. The basic assumptions of the approach remain the same, however the algorithm 
has to be modified in some essential aspects. First of all, one has to build the mechanism of temporal 
derivative approximation into the numerical procedure together with necessary measures to stabilize the 
algorithm. We will propose two methods of doing this, constructing in fact two different solvers. The 
second fundamental difference between the self-similar and time-dependent formulations is that in the 
latter case, the crack length L{t) becomes now an element of the solution, which should be looked for 
simultaneously with the crack opening w(i, x) and the first term of its asymptotics near the crack tip 

Wo{t). 

The basic system of equations for the transient problem is composed of: the governing equation (|22p . 



11 



the boundary conditions (jl6p . the initial conditions (|17p and the integral equation defining the crack 
length (mi). 

To avoid using multiple subscripts let us adopt the following manner of notation: 

w{tj,x) = w{x), 'w{tj+i,x) = W{x). (45) 
Consequently, the asymptotic representations of the solution near the crack tip read: 

w{x) = wo(l — xY/^ + Aw, 

a; ^ 1, (46) 

W{x) = Wq{1 - xY'^ + AVF, 

where C, is defined in ([55)1 and 

/^w^O{{l-xf), = 0{{l-xf), as x^l. 

3.1 Solver using self-similar algorithm ( I40p — (14 Ih 

Below, we show the way to convert the initial boundary value problem defined by equations (|16[) . (j2ip . 
(|22|) to the form which may be tackled by the integral solver in the form (|4T|) for the self-similar solution. 
Obviously, system (|1T|) is to be supplemented with an additional equation defining the crack length. 

The main idea of the approach is to use the temporal derivative as one of the dependent variables in 
the numerical procedure. To achieve this, we compute the derivative of the solution at each time step 
t = in the following iterative process: 

dt V dt . ^^^^ 

At ^\ ) dt ' 

where superscripts refer to the number of iteration, At = tj+i ~ tj and the values of cr'-'"'"^^ are to be 
defined later. The first approximation of the temporal derivative is 

Note that the derivative at initial time t = can be immediately found from the governing equation 
([^^ by substitution of the initial conditions. 

Substituting (|47)) into ((22|) one obtains the problem ((26|) and p7| with respect to the unknown 
function W{x) by exploiting the following convention: 



^ At 

Ql = -W + 



('+1) 



q (i+l) , ,2 



At 
and 

q*it) = 3qo{t)W^'/^L(^l (50) 
Finally, the crack length from equation (PT|) can be iterated as: 



(51) 
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Here the integral from ([?!]) is computed by the trapezoidal rule, taking into account the information 
obtained from the previous time steps (t < tj). 

In this way the governing partial differential equation and all the boundary conditions have been 
transformed to the formulation used previously in the self-similar problem. As a result, respective 
integral operators ([M)) - p7p and ([5^ remain the same with modified arguments ([^^ and ([50)1 . where 
([5T|) should be taken into account. 

To choose the value of the parameter (t(*+^) at each iterative step, let us recall that the best perfor- 
mance of the algorithm described in subsection 12.21 has been achieved near /3 = 1/3. Taking this fact 
into account in (PU|) 9. one can choose 



At 



9(L(*+i)) 



(52) 



Thus, the inverse operator for A, defining the solution of the transient problem at the time step 
t = tj+i, can be expressed in the following manner (compare (j40p ): 



Wo, L, AW, 



dW 



A-'{l/3,q:,q*o) 



(53) 



The iterative algorithm of the solver can be described by the system: 

'G3 {l/3,W^'+'\AW^'\q*o) =0, 
AVl^('+i) = Gi (l/3,Wo(*+'\Aiy«) + 
G2 (l/3,W(j''+'',Al^W,q; 



(54) 



L(*+i) = Gg ( W, 

dW^'+^^ 
~dt 



G4 



dW^^^ 
dt 



Note that parameters = q^'"'' , ql ~ Q;*'"'^^'' and (t('+^) are also iterated according to (051)1, (|50p and 
([5^ . To finalize the algorithm, it is enough to define the initial forms of the crack opening and 
the crack length L^^\ This is achieved by choosing 



(55) 



which finishes the description of the computation scheme for the fixed time step At = tj^i — tj. Let 
us recall that the temporal derivative of solution at time t = is computed using the initial conditions 
(fT7)) . while for any next t = tj we utilize representation (|Tf)) . 

We would like to underline here the fact that as the output of the proposed algorithm, one obtains 
not only the solution of the transient problem, that is, the crack length L(t) and the crack opening 
w{t,x), but also the temporal derivative of the latter, w'^ltjx) and the crack propagation speed, Vo(t), 
from (1201). 



3.2 Solver based on a modified self-similar algorithm 

Temporal derivative of the solution can be taken in the following form: 



dW 
~dt 



W 



dw 

~dt' 



(56) 
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which gives the error of approximation of the order O(At^). Note that any other two-points finite 
difference definition yields only 0(At). 

Unfortunately, a direct use of the algorithm formulated in section 12.11 is not, generally speaking, 
possible, as it may fail for small time steps. Indeed if one takes tr = 2 and sufficiently small value of At 
in the representation (|49p 9. then the value of the parameter /3 may be far away from its optimal magnitude 
/3 = 1/3. This in turn, would at least result in a deterioration of the efficiency of computations. 

On the other hand, it is quite clear that boundary value problem pG]). ((27| is solvable for large 
values of the parameter /3. Indeed, performing asymptotic analysis, one shows that the solution can be 
represented in the form: 

u{x) ^ ~q'l{x) +boix) +bi{x), a; 6 (0,1), (57) 

where bo and &i refer to two boundary layers, accounting for the boundary conditions (P7)) i and ([?7)) 9. 
respectively. 

Thus the algorithm should be modified to be able to deal with large values of (3. To achieve this goal, 
let us look at the original representation ([M)) . where only the last term in the right-hand side depends 
on /3. This term violates the convergence of the iterative process for large /3. To prevent this from 
happening, at each iteration we supplement the term in question with an auxiliary so-called 'viscous' 
term in the form V{x) = /3uq(1 — x){Co + Ci(l — x)), where the constants, Co,Ci, are computed by 
comparing both the original and the viscous terms. In this way one can construct a modified algorithm, 
schematically represented in the following manner 

0, 



(58) 



= G 



5 [Wo 



(' + 1) 



At 



where 



Qi = 



-w 



At 
~2 



dw 

"dt' 

dw 

2 



3 



(59) 



The function qo{t) is defined in the same way as previously (see (jSOp V At each time step the initial 
value of W is taken in the form (|55|) . We do not show here the precise definition of the operator G12, as 
it may be easily derived by merging operators Gi and G2 with the viscous term V{x), where respective 
constants are computed, for example by the least squares method. 



3.3 Analysis of the algorithms performance. 

The aim of this subsection is to analyze and compare the performances of the algorithms formulated 
in subsections 13.11 and 13.21 To this end, the benchmark solution used previously for the self-similar 
formulation is utilized, for the time dependent term %lj{t) = (1 + t)'' (compare ([55)) ). First tests are 
performed for 7 = 1/5. 

In the following, the notations solver 1 and solver 2 are attributed to the algorithms ([5^ and ([55]) . 
respectively. 

Let us analyze the influence of the spatial mesh density on the accuracy of computations done by 
both solvers. For this reason, we start with a single time step (in our case At = 10~^) and carry out the 
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computations for different numbers of nodal points TV, ranging from 10 to 100. The following parameters 
are used for the comparison: the maximal relative error of the crack opening, Sw, the relative error of 
the crack length, SL, and the maximal relative error of the temporal derivative of the crack opening, 
Swt- 

The results of computations are depicted in Fig. HJi). It shows that regardless of the considered 
parameter, solver 2 always provides better accuracy (two orders of magnitude for the crack opening, w, 
and an order for the crack length, L). The only exception is for Swt, whose values are of the same order 
and solver 1 may even give a bit lower errors. However, this does not result in a better accuracy of the 
two remaining components of the solution: w and L. Note that for the first time step, the accuracy 
of computation of the crack length, L{t), is of two orders of magnitude better than that for the crack 
opening, 5w, regardless of the solver type and the number of the nodal points, N . 

Similarly to the trend observed in the self-similar formulation, computational errors stabilize for some 
critical value N = N^.{At). Surprisingly, for the transient regime where the temporal derivative plays a 
crucial role, the critical iV* appears to be even slightly lower than that for similarity solution. Thus, for 
both solvers it is sufficient to take only 40 points to achieve the maximal level of accuracy. 




10" 
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-8L 
Swt 
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Figure 4: The errors of computations: 6L,dw,dwt as functions of a) spatial mesh density (number of 
nodal points, N) for the fixed time step At = 10^^; b) time step At for fixed number of the nodal points 
= 40. Line with markers refer to solver 2 



Remark 2. From first glance it may be surprising that the accuracy of the temporal derivative, 5wt, 
is up to four orders of magnitude worse than that for bw in the case of solver 2. However, this fact can 
be easily explained when one analyses the estimation: Swt ~ 2wSw / {Atwt) which follows from (|56p for 
small At. Computing the multiplier in the right hand side of the estimation, one obtains the respective 
value of the order of 10"'. Note that the recalled formula docs not account for the error of the method 
of Wt approximation itself. In some cases (very small time steps) this value may be comparable to the 
former, and thus essentially infiuence the overall error. 

Remark 3. At the first time step, there also exists a direct relationship between the error of the 
crack length, SL, and the multiplier of the leading term of the asymptotics (^S)) . SWq. For small values 
of At it reads: 26L = W^SWoAt. Having this relation we do not show a separate analysis for Wq and 
consequently, Vq. 

All the results presented above were obtained for a single time step At = 10~^. In order to illustrate 
the influence of At on the solution accuracy, a number of computations were done in the interval 
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At € [lO""^,!] for a fixed number of nodal points, N — 40. The results describing the solution errors 
are shown in Fig. |4]b). One can conclude from them that the errors decrease when reducing the time 
step (at least in the analyzed range of parameters). Simultaneously, one can expect that this tendency 
should have its own limitation for a fixed number of the nodal points, N, and starting from some small 
value At{N), it will reverse to the opposite. For both solvers, a fast decrease of SL is observed as At 
gets smaller {SL « alO~*At^ and a ^ 1). 

Now let us analyze the accuracy of the solution on a given time interval t G [0,ix]- As, it follows 
from the results presented in |17| and [22j . after the initial growth to the maximal value, the error of 
computations stabilizes at some level or even decreases, and thus further extending of the time interval 
does not contribute to the deterioration of accuracy. Bearing this in mind we take here Ik ~ 100. First 
we test how the number of time steps K affects the accuracy of computations within the same time 
stepping strategy. For both solvers we took a fixed number of spatial mesh points, N ~ 40, changing 
the number of time steps K from 10 to 300. The utilized time stepping strategy was the same as that 
used in [52] for the Crank-Nicolson scheme. It is given by the following equation (i = 1,2,..., Ky. 

^,:^(^-l)-5t+^i^^^i^^(^-l)^ (60) 

where bt is a parameter controlling the first time step. Note that by increasing the value of K one 
distributes the time points near i = almost uniformly. 




10 50 100 150 200 250 300 

K 



Figure 5: Influence of the number of time steps, within time stepping strategy on accuracy 

of the computations for the fixed non-uniform spatial mesh (N = 40, = 3). Lines without markers 
correspond to solver 1, lines with markers refer to solver 2. 

The achieved accuracy, as a function of the number of time steps if, is presented in Fig. [SJ It shows 
that for the same value of K solver 2 provides much more accurate results again and this advantage 
increases with growing K. One can expect, as it follows from the discussions after Fig. |4] which refers to 
a single time step, that each solver gives the solution errors satisfying the estimation: SL << Sw << Swt 
for any set of the input parameters. Indeed, solver 2 supports this statement as can be seen Fig. [SJ 
while for solver 1 surprisingly another trend is observed: 5w < 5L « dwt- 

To explain this apparent paradox, let us recall that all errors in Fig. [S] are taken as the maximal 
values over the time-space domain. However, since for the assumed time stepping strategy the time 
steps increase with growing time, one can expect (compare the results in the Fig. HJd)) that at some 
point 6L may become greater than Sw. Moreover, the error accumulation in successive time steps may 
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Tabic 1: Comparison of the results obtained for the solver developed in [35] and two integral solvers: 
solver 1 and solver 2. The same benchmark solution with the leak-off vanishing near the crack tip is 
considered. 



additionally influence the relationship between 5L and 5w, especially if they are close to each other as it 
is in the case of solver 1. Fig. [5] depicts the evolution of computational errors in time. For both solvers 
dw and 5wt reach their maximal values and stabilize or decrease for t < tx- In the case of solver 1 one 
can observe an intensive error accumulation for the crack length, 5L. Indeed, it can be seen from the 
Fig. [51 that there exists a moment when the error curves for 5L and 5w intersect for the solver 1. In 
the case of solver 2, respective errors of solution ((5i, Sw and 5wt) have essentially different values. As 
a result, this effect does not take place in the time interval under consideration. However, it may be 
encountered for Ik > 100, as suggests the trend observed for t close to Ik- 




t 

Figure 6: Distribution of the solution errors in time. Line with markers refer to solver 2. The time 
stepping strategy ([60]) is used for K = 30. The number of nodal points, TV = 40. 

It is interesting to compare the solution errors generated by the two integral solvers with those 
obtained in [^^ for a dynamic system (DS) approach. In the latter case a standard MATLAB solver, 
odelSs, was employed. The best results were obtained for a uniform spatial mesh, the presented data 
corresponded to = 100. The time stepping strategy chosen automatically by the solver (whose 
character is approximated by ((60|)) accepted 242 time steps. The maximal relative error of the crack 
opening, Sw, and the crack length, 6L, were 5.0 - 10"3 and 6.8 - 10-3 (see Table [T]), respectively. 

When analyzing the data collected in Table [TJ one concludes that it is sufRcient for any of the 
integral solvers to take only A^ = 40 nodal points and if = 30 time steps, as suggested the previous 
analysis, to have better (but comparable - solver 1), or much better {solver 2) results. Indeed, the 
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corresponding maximal errors for the integral solvers are: Sw ~ 3.7 • 10~^, SL ~ 5.2 • 10"'^ for solver 
1, and Sw = 1.1 • 10~^, SL = 2.4 • 10~^ for solver 2. In other words, the first solver provides the same 
accuracy for the crack opening and the crack length as the DS solver using much greater numbers of 
nodal points and time steps, while the second one, under the same conditions, improves the results at 
least one order of magnitude. It also shows that the solver 2 yields one order of magnitude better 
accuracy of the crack propagation speed, Vb, than the solver 1. 

The new algorithms allow us to automatically compute the temporal derivatives in the solution 
process. The respective errors, Swt^ are: 7.5 • 10~^- solver 1 and 1.5 • 10~^ - solver 2. We decided 
to compare these figures, with the ones obtained in postprocessing (here, also the DS approach was 
examined). To this end two FD schemes (2-points and 3-points) were used. This time, the corresponding 
errors, Swt, were: 4.6 • 10~^ and 2.0 • 10~^ for solver 1, 4.6 ■ 10~^ and 3.4 • 10~^ for solver 2 and the same 
value 4.1 • 10~^ for both schemes in the case of DS. It is worth mentioning that the values obtained for 
DS appeared at the first time step. Then, the errors decreased with time and stabilized to give their 
minimal levels of 10~^ and 4.8 • 10~^, correspondingly. 

As can be seen, the integral solvers give at least one order of magnitude better accuracy of wt than 
the DS. Moreover, while the postprocessing gives smaller (but comparable) error for the solver 1, solver 
2 returns more accurate values of Wt than those obtained in the postprocessing, even for the 3-points 
FD. Finally, apart from the fact that Sw and SL for the DS and solver 1 look comparable in values, the 
quality of the computation is better for the new solver as is clear from the postprocessing analysis. 

Just for comparison we also present in Table [T] the results obtained for the spatial mesh composed of 
only five nodal points, N ~ 5. It turned out that even for such a drastic reduction of the mesh density, 
the solution accuracy for most of the parameters is of the same order as for N = 40. Interestingly, solver 
1 exhibits almost no sensitivity to this mesh reduction. In fact the distinguishable differences can be 
observed only for the solver 2. 

Let us now analyze the distribution of crack opening error, Sw, in time and space. The respective 
results are presented in Fig. [7] a) - solver 1, and Fig. [7]b) - solver 2. In both cases the maximal errors 
are located at the crack tip while the error distribution in time follows the trend visible in Fig. [SI 




Figure 7: The relative error of the crack opening obtained for = 40 (nonuniform mesh, g = 3) and 
the time stepping strategy ((60)) with K = 30 . Fig. [7^) corresponds to the solver 1, while Fig. [7}d) refers 
to the solver 2. 

Fig. |5] shows the distributions of Swt obtained by the integral solvers. It confirms our previous 
observation, that solver 2 always provides better results than solver 1. Moreover, the greatest error in 
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the case of solver 1 is located at the crack inlet and the lowest at the crack tip, while solver 2 gives 
approximately the same values of 5wt along the crack length. 




Figure 8: The relative error of the temporal derivative of the crack opening. Solution obtained by: a) 
solver 1, b) solver 2 for N ~ 40, non-uniform mesh {g = 3) and the time step strategy (pO)) with K ~ 30. 
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Figure 9: The errors of solution for two variants of 7: Fig. [S^) 7 = 0, Fig. I^b) 7 = 1/3. Line with 
markers correspond to solver 2. 

In the last test in this subsection we analyze the relation between the regimes of crack propagation 
and the performances of respective solvers. As mentioned previously, the benchmark solution in form 
(|29p can be used to imitate various dynamic modes of the crack evolution. So far we have utilized the 
exponent of the time dependent term of the value 7 = 1/5 which refers to the constant injection flux 
rate. Now, let us consider two other variants of 7: i) 7 = - for this choice the normalized crack opening 
is constant in time; ii) 7 = 1/3- this value corresponds to the steady state propagation of the fracture. 
For the computations the same spatial mesh as before is taken {N = 40). The number of time steps 
accepted within strategy ([60|) . K, ranges from 10 to 300. In this way the graphs (Fig. [9li)-b)) describing 
the solution errors in the function of K were prepared for both values of 7 (similarly as in Fig. [5]). The 
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analyzed accuracy parameters were: relative error of the crack opening, Sw^ relative error of crack the 
length, 5L, relative error of the crack opening temporal derivative, Swt^ for 7 = 1/3 and absolute error 
of the crack opening temporal derivative, Awt, for 7 = 0. 

The results depicted in Fig. ^ show that for 7 = one obtains much more accurate results than for 
7 = 1/3 that could have been predicted {wt = 0). However, it is a surprise that for 7 = solver 1 
provides better solution accuracy than solver 2. Although the difference is moderate in case of Sw and 
SL, the values of Awt vary by at least two orders of magnitude. From Fig. [9^) it follows that for this 
regime of crack propagation the solution accuracy cannot be improved by simple refining the temporal 
mesh, and for solver 2 even a reverse relation is observed. 

The situation is quite different for 7 = 1/3 (Fig.[3]D) ). This time again solver 2 proves its advantage 
over solver 1 for all the analyzed parameters. For solver 2, it is sufficient to take only 30 time steps to 
have much better results than those provided by solver 1 for 300 steps. The solution accuracy can be 
improved by increasing the number of time steps, however it seems that for solver 1 the saturation level 
is close to if = 300. A simfiar trend was observed for 7 = 1/5 (see Fig. [5]). 

A direct conclusion from this test is that for different modes of crack propagation there are different 
optimal time stepping strategies. This should properly accounted for especially in the cases when the 
values of injection flux rate or leak-off to formation change appreciably in the considered time interval. 

The aforementioned analysis prove that in terms of accuracy in most of the cases solver 2 is much 
better that solver 1 with respect to all computed components of the solution: the crack length, L, 
the crack opening, w, its temporal derivative, Wt and the fracture propagation speed, Vq. However, 
for some regimes of crack propagation (low values of 7) solver 1 may give comparable or even slightly 
better results than solver 2. The advantage of solver 1 is better efficiency of computations: the time of 
computations for this solver was on average one third lower that for solver 2. 



3.3.1 Example with singular leak-off regime. 

In ([3]) we have assumed that the behaviour of the leak-off function near the crack tip can be described by 
a power law, giving in the worst case a square root singularity. Such a limiting behaviour corresponds to 
the Carter leak-off model [5]. As a result, although the leading term of the asymptotic expansion for the 
crack opening near the fracture tip remains the same, the higher terms change, disturbing the solution 
smoothness (see [15| ) . A comprehensive analysis of this case was done in jl7j where it was also proved 
that the deterioration of solution accuracy can be prevented by employing the second asymptotic term 
in the computational algorithm. 

In this subsection we show that the algorithms developed in the paper are capable of tackling this 
kind of problems without any additional modifications. To this end let us consider another benchmark 
solution u{t,x) (see (|63| ) defined by the functions: 

h{x) = {l-xY^^l + s{x)), s(x) = i(l-x)i/6, (61) 

5 

with the same function ipit) with 7 = 1/5. 

One can easily check, that the above form of s{x) results in a singular behaviour of qi, with the 
leading term of the order 0((1 — x)^^/^) as x — > 1. The value of multiplier uq in ((63)) was taken in such 
a way to make the benchmark comparable with the one used previously, in a sense of an average particle 
velocity. Indeed, in [52] for a fluid velocity defined as V = q/w, a parameter describing its variation 
along the crack length was introduced: 



lv{t) = [max{V{x,t)) — ■min{V{x,t)) 



V{x^ t)dx 



(62) 



This parameter refiects indirectly the balance between the fiux injection rate and leak-off to formation. 
It was also shown there that it has a decisive influence on the accuracy of computations (the greater 
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value of 7t,, the greater error of the computations). For the benchmarks with comparable values of 7^, 
one can expect similar accuracy of the computations. This trend was also confirmed in |17| . In our case, 
the deterioration of the solution smoothness near the crack tip is an additional factor which contributes 
to the increase of the computational error. 

Note that because of the chosen structure of the value of 7t, is constant in time for all benchmarks 
considered in this paper. For the benchmark P5)) 7^ yields 0.408, while (|5T|) one leads to 7t, — 0.411. 

The computations were done for the same spatial (A'' = 40, g = 3) and temporal {K = 30, strategy 
([50]) ) meshes as previously considered. The distributions of the errors 5w and Swt are shown in Fig. [TU] 
- Fig. [Til respectively. Comparing these results with those obtained for the finite leak-off regime (see 
Fig. [7] - Fig. IS]) one can see that the accuracy of computations decreased significantly. 




Figure 10: The relative error of the crack opening. Solution obtained by: a) solver 1, b) solver 2. Other 
parameters in the computations were: N = 40, K = 30, g — 

In the case of solver 1, the crack opening error, Sw, shown in Fig. [10] a) is of one order of magnitude 
greater than that reported in Fig. [7] a). For solver 2 in turn, 5w depicted in Fig. [10] b) is almost two 
orders higher than the one obtained previously for qi vanishing at a; = 1 (Fig. [7]b)). A pronounced 
jump of the error is observed at the crack tip, especially for solver 2, which explains the deterioration 
of accuracy when comparing with the finite leak-off case. However, if one considers the accuracy of the 
solution away from the crack tip, it is still of the same quality as before. 

The same trend can be observed for the temporal derivative of the crack opening, Wt, see Fig. [11] 
Opposite to the benchmark case (|43p with the finite leak-off. the distribution of Swt for solver 2 becomes 
highly non-uniform, with distinct increase at the crack tip. 

This deterioration of the solution accuracy near x — 1 for the Carter leak-off model should not be a 
surprise, as the algorithm described above in p4)) - (|37]) accounted directly only for the first (leading) 
term of the asymptotic expansion for the crack opening. For the singular leak-off model, the function 
Am and other integrands employed in Gi and G2 are not sufficiently smooth near the crack tip, which 
increases the errors of integration. However, one can counteract this tendency by accounting for further 
asymptotic terms in the algorithm, which will be illustrated in the following with the example of two 
terms of w expansion. 

Let us investigate the evolution of solution errors {6w, SL and 6wt) in time in the way it was done in 
Fig. [6] for the case of the finite leak-off. In Fig. [12] a) we show the results obtained by direct execution of 
the algorithm ([34]) - ([37]). Fig. [12] b) depicts the case of a modified algorithm employing also the second 
asymptotic term of the crack opening. When analyzing the results, one observes again the trend of SL 
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accumulation for both solvers. Interestingly, SL has the same level for both the original and the modified 
algorithms. In the case of solver 1 for other analyzed parameters as well only a very slight improvement 
is reported. However, for solver 2 one obtains a pronounced reduction of computational errors for w and 
Wt- To have a further advance here, the next terms of asymptotic expansion of the crack opening may 
be taken into account. On the other hand, for the Carter leak-off model the number of nodal points for 
which the accuracy saturation is observed is larger than for a finite leak-off. Thus, if necessary, one can 
also improve the accuracy by increasing the value of N and by taking a greater number of time steps K. 




Figure 11: The relative error of the temporal derivative of crack opening. Solution obtained by: a.) solver 
1, b) solver 2. Other parameters in the computations are: = 40, if = 30, g = 3. 




t t 

Figure 12: The relative error of the temporal derivative of crack opening. Solution obtained by: a) 
original algorithms, b) modified algorithms. The lines with markers refer to solver 2 

From the analysis given in this subsection two main conclusions can be drawn. First, that the 
problems with singular leak-off regimes can be directly tackled by the proposed algorithms. Also in 
such cases solver 2 yields more accurate results than solver 1. The second conclusion is that in order to 
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have the solution accuracy comparable to that achieved for non-singular Icak-ofF, one has to use a larger 
number of nodal points and/or employ further terms of asymptotic expansion for w in the algorithm. 

4 Conclusions 

Wc would like to itemize the following conclusions as a resume of this paper: 

• Presented approach can be efficiently used for tackling the PKN model of hydrofracturing and may 
be adopted for multifracture systems. 

• Both new solvers provide better computational accuracy than the conventional algorithms from 
P^ . Moreover, comparable accuracies can be achieved here at much lower computational cost, as 
the new solvers enable us to drastically reduce the densities of spatial and temporal meshes. 

• New solvers are appropriate for directly tackling the problems with different fluid flow regimes, 
including various injection flux rates and singular leak-off. 

• In order to increase the efficiency and accuracy of computations, it is advisable to employ at least 
two asymptotic terms of the crack opening, w. 

• The developed algorithms do not require any regularization techniques. The boundary conditions 
are imposed directly into the numerical scheme. The speed equation plays a crucial role in the 
analysis. 
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Appendix A: Numerical benchmarks 

Let us define a set of benchmark solutions useful for testing different numerical solvers. Consider a class 
of positive functions C-|-(0, 1) described in the following manner: 

C+(0,1) = {h e (7^(0,1) nC[0,l], 

lim (1 - xy^^'^hix) = 1, h{x) >0, x£ [0, 1)}. 

By taking an arbitrary h G C+(0, 1), one can build a benchmark solution for the normalized formulation 
of the problem as: 

u{x) — UQiljj{t)h{x). (63) 
where functions i/jj{t) and h{x) are specified below. On substitution of (|63| into (|26|) one finds: 

qi{t,x)= (64) 
juo[^ (^xh''ix)h'{x) + 3{h^x)h'{x)y^ -h{x)]i^^{t), 

where two sets of the benchmark solutions can be considered. For the first one we choose V'i(0 ~ 
and /3 = 2/3, a = 1, while for the second, V2(i) = (a + t)'^ and jS ^ 27/(37 + 1), a = (37 - l)/7. 
Corresponding crack lengths are defined in (|25p and pop . respectively. Finally, the injection flux rate is 
computed from the boundary condition (|16P i : 



Qoit) - -^j^h^O)h'{0), (65) 



while the initial condition reads W{x) — uo'iJjj{0)h{x) 
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Note that when taking the function h{x) from the class C+(0, 1) in the following representation: 



s G C2[0, 1], s{x) > -1, xe [0,1) 

qi automatically satisfies the condition qi{t, x) = 0((1 — a;)^/^), x — ?► 1. 

The presented benchmarks allow one to test numerical schemes in various fracture propagation 
regimes (accelerating/decelerating cracks) by choosing proper values of the parameter 7 (see [12]). 
Additionally, if one reduces the requirements for the smoothness of s{x) near x = 1, assuming that 
h e C^[0, 1) Pi iJ"(0, 1) the benchmark can serve to model singular leak-off regimes (compare with (pT|) ). 

Note that the zero leak-off case cannot be described by the aforementioned group of benchmarks. 
However, an analytical benchmark for this regime, represented in terms of a rapidly converging series, 
has been developed in [3T]. 
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